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1 . 


SUMMARY 


A research effort was conducted with the goal of reducing 
computer time of a Navier Stokes Computer Code for prediction of 
viscous flow fields about lifting bodies. A two-dimensional, time- 
dependent, laminar, tran&onic computer code (STOKES) was modified to 
incorporate a non-uniform timestep procedure. The non-uniform time- 
step requires updating of a zone only as often as required by its own 
stability criteria or that of its immediate neighbors. In the 
uniform timestep scheme each zone is updated as often as required by 
the least stadale zone of the finite difference mesh. Because of less 
frequent update of program variables it was expected that the non- 
uniform timestep would result in a reduction of execution time by a 
factor of five to ten. Available funding was exhausted prior to 
successful demonstration of the benefits to be derived from the 
non-\iniform timestep method. 
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2.0 


INTRODUCTION 


Two elements lie at the core of any satisfactory numerical 
procedure for computing laminar and/or turbulent transonic airfoil 
flows, namely, (a) a general method for integrating the equations of 
continuum mechanics, to which equations governing the evolution of 
turbulence variables can be added, and (b) a mathematical model that 
describes the fluctuating component of an arbitrary turbulent field. 

To the extent possible, both the numerical method and the turbulence 
model should offer proven reliabi lity; the choices proposed are believed 
to best meet such a criterion, r.amely: 

1. A flow field computer code, which integrates 
the Navier-Stokes equations in two spatial 
dimensions, and 

2. The Saffman Model of turbulent flow. 

For each of the many laminar flows (References 1-7) to which it has 
been applied, the presently operational flow field computer code has been 
shown to solutions accurate enough for ainy engineering purpose . 

Furthermore, in the two-dimensional test calculations completed to date, 
the less proven Saffmein turbulence (Ref. 8-9) model, in combination with the 
flow field computer code, has provided cin accurate solution to the turbulent 
flat plate problem, (Ref. 10) a turbulent boundary layer-shock wave inter- 
action (Ref. 10) and tne separation phenomena that occur on a compression 
corner in supersonic flow. (Ref. 10) 

A version of tlie Afton plane '-ymmetric code (A2P) was developed 

for the NASA Ames Research Center to solve a laminar transonic airflow 

problem. The flow field was computed about the lifting NACA 64A410 air- 

G 

foil at Mach .72, Reynolds number (based on chord length) of 1.75 x 10 

and at an angle of attack of 4 ® (ref 11). Starting from the impulsive conditions the 

calculation was carried out to the time it takes a freestream particle 
to travel one chord length. The calculation required slightly less thein 
five hours of CDC 7600 computer time. In the course of tlie calculation 
the flow separated on the leeward trailing edge, a leeward compression 
wave formed, and approximately 50% of the measured lift was computed. 

( 
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However, it was estimated that an additional ten hours on the CDC 7600 
computer woald be required to complete the problem. 

Conputing times of the order of fifteen hours on the CDC 7600 
to complete one transonic airfoil calculation make transonic flo>'J 
computations about lifting bodies impractical for design purposes. A 
speed-up factor of the order of 10 would be necessary for numerical 
computation of treuisonic flow to be in the practical range. Coupled 
with a better estimate of the initial conditions than impulsive inital 
conditions, and with a speed-up factor of ten, transonic lifting airfoil 
problems cam be solved in less theui an hour on the CDC 7600 computer. 

The main objective of the work conducted in this program was 
to develop a computer code wich would be approximately an order of 
magnitude faster tlian the present computer program. This was to be 
accomplished by use of a "Non-Uniform Timestep" in the computation 
procedure to improve computation efficiency. 
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3.0 SYMBOLS 

E Specifie internal energy 

j Index specifying streaunline-like direction 

k Index specifying potential -like direction 

M Mass 

n Time index for finite difference equation 

p Pressure 

At Tiroes tep increment 
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4.0 NON-UNIFORM TIMESTEP FOR IMPROVMENT IN COI'IPUTATIONAL EFFICIENCY 

OF NAVIER STOKES SOLUTIONS 

The STOKES computer code was developed for the NASA Ames Research 
Center from eui AFTON pleme symmetric code (A2P ) . It is a two-dimensional 
time-dependent code used to predict laminar transonic flows eibout lifting 
bodies. In common with other AFTCW codes, a uniform timestep was incor- 
porated in the original version of the program. The uniform timestep 
procedure requires that all of the variables of motion in every zone of 
the mesh be updated or calculated as often as required by stability 
considerations for the least stable zone of the mesh. The timestep is 
dependent upon the time required for the speed of sound to cross a zone. 

It was anticipated that significant computer time savings wf uld be achieved 
if the non-uniform timestep procedure were incorporated in the code. This 
allowed the varicd)les of motion to be calculated for a give zone of the 
finite difference mesh only as often as required by considerations of 
stability of that zone and its immediate neighbors. It was estimated that 
a speed-up factor between five and ten would result from the incorporation 
of non-uniform timestep procedure. 

4.1 Treatment of Timestep Intermesh Boundaries 

One of the most crucial aspects of the non-uniform timestep procedure 
relates to the way boundaries between adjacent zones having different times teps 
are treated. In order to maintain AFTON conservation properties across the 
interface it had been previously determined by Trulio (Ref. 12) that correct time 
centering of the pressures used to update the momenta was essential. As initially 
conceived, the interface treatment would assume that the adjoining large /.t zone 
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would be treated in nearly all respects as the small timestep level zone 
(a micro zone) amd updated with a small timestep. In this 
macro boundary zone the explicitly calculated pressure available from 
the first iricrostep update (there would be two such calculations in order 
to bring the zone correctly forward in time equivalent to a single 
macrostep calculation) is the correct pressure to use in updating the 
adjacent momentum zones. This pressure would be used for both tr.lcrostep 
momentum updates amd the single macrostep momentum update. This procedure 
would have correctly treated the momentum conservation constraints in the 
tramsition from a micro zone to a macro zone. 

The method described above does not correctly conserve the total 
energy in the same tramsition region. For a macro level zone, j, in the 
Lagramgii\n case, 


M. 





( 1 ) 


must be true toconser^ total energy; where is the zone mass, the 
specific internal energy, the zone pressure amd 3^^ the chamge in 
volume. Superscripts refer to centering in time; n is a microstep 
increment. Macrozone energy is updated in time from n - 1 to n + 1, an 
interval of 2n or two microsteps. Had we treated this boundary macro 
zone as a micro zone, as described in Paragraph 2, two such updates would 
be required. 


M . (E . - E . ) * p . AV . 

1 3 D 3 3 

M. (e""^^ - e’’ - -p" Av” 
3 3 3 3 3 


( 2 ) 

(3) 


The pressure would no longer be correctly centered in relation to the 
internal energy. The resulting change in internal energy from E^ ^ 
to E would not be the same as that calculated by Equation (1) . 

The interface was revised to incorporate a different method which 
preserves the general conservation properties. Macro zones in the interfaces 
were to be treated as macro zones in general respects; em intermediate 
pressure centered at n is calculated in the macro zone and used to update 
the adjoining momentum zones. This pressure at n is used with Equation 
(1) to explicitly calculate the internal energy at n + 1 and so preserve 
total energy conservation. Logic and coding were developed and inserted 
in the updated program. 
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4.2 Development of Non-Uniform Timestep Logic 

The logic for controlling the zone propertiesupdating sequence 
and the assignment of timestep levels for thermodynamic and momentum 
zones was included in two new subroutines inserted in the code. These 
subroutines, ZONEL and TMESTP are shown in the program flow diagram 
presented in Figure 1. ZONEL is the routine which assigns timestep 
levels to the thermodynaunic and momentum zones. It is called from 
FLOW at the beginning of a new macrocycle. At that point all zones 
have been updated to the same time and new levels must be assigned for 
the new updating sequence. 

The TMESTEP subroutive calculates and stores the actual timestep, 
based upon sted^ility criteria,- for each zone. It replaces coding formerly 
in MVS. This routine is also called from FLOW, but after MVS and after 
the final macrolevel update. The timesteps are then available for 
ZONEL to perform its required testing and assignment at the beginning 
of the next macrocycle. A listing of subroutine ZONEL is included as em 
Appendix to this report. 

The progr^t..l, as written, is currently restricted to only two 
timestep levels, although additional t: e reduction should result from 
multiple time levels. Thus, zones which are at least twice as i.aige as 
the smallest zone require updating only every other cycle. Zones which 
are four times larger than the smallest zone. would normally require 
updating only every fourth cycle, etc. A temporary restriction in 
timestep levels was inserted to simplify checkout emd debugging. 

Coding is in place which will eventually allow up to five timestep 
levels . 

4 . 3 Development of Conservative Logic 

In order to maintain the conservative aspects of the AFTON 
numerical method an approach was tedcen which utilized momentum and 
mass accummulators to determine mesh boundary fluxes. This approach, 
because of its inherent simplicity, offered s'jbstemt ial savings in 
execution time in addition to the savings which were expected to result 
from the non-uniform timestep. 
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Figxire 1: General Flow Diagram for Computer Program With Non-Uniform Timestep 

Scheme. 
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Conplete conservation also required calculation of the work 
done on the surfaces of any body placed in the flow field. The AFTON 
differencing scheme results in a half momentum zone around each of the 
mesh boundaries and around any body in the mesh. Logic had to be im- 
plemented which accounted for the fact that the work calculation for 
the body surfaces included the momentum trcinsported across these half- 
zone boundaries. 

The various difficulties which were encountered while incor- 
porating the conservation logic eventually led to exhaustion of the 
available funding for this program. It is estimated that about one- 
half of the conservation logic had been inserted in the code. Sub- 
stantial additional programming and debugging effort is required to 
con^lete thi s task and insure that all of the conservation checks 
of the AFTJ?' T.echod are preserved. 

4.4 Checkout and Debugging of Non-Uniform Timestep Coding 

A test problem was run tc check out the code for leuninar 
flow over a right circular cylinder. Three cycles were run with time- 
step restricted to two timestep levels. The program appeared to execute 
normally. However, validation of the changes requires completion of the 
conservation logic emd coding. 
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5.0 


CONCLUDING REMARKS 


An effort weis undertaken to incorporate a non-uniform times tep in 
a two-dimensional, time -dependent, Navier-Stokes solution of the laminar, 
transonic flow over a lifting body. Available funding was exhausted 
during the implementation of the conservation logic. Without the conservation 
properties the accuracy of the numerical calculation is uiicertain. Successful 
demonstration of the non-uniform timestep as a method for improving the 
efficiency of Navier-Stokes solutions requires additional programming and 
debugging effort. 
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270 


Jl ■ JB 

20NEL.99 

2 ; 1 


J2 X JB « 2 

20NEt.90 

27i 


JJ * JB 

ZOnEl.JI 

27« 

641 

continue 

2JNEL.92 

27« 


JNL(JI.L) > »INO(J2L(JJ,L).J2L(Jl.U5}>J2u(JZ>un 

/jNEw.PJ 

Jl J 


IECJI.nE.JB) CQ to 66 

23NEL.94 

SIS 


Jl» JPn 

2Jnel.9> 

S16 


J2< JBn t 1 

ZJNEU.96 

J17 


JJ» JB« • 1 

2JNEU.67 

llO 


CO to 641 

ijNEL.99 

iiO 

65 

continue 

20NEL.96 

liO 


Jls JB 

2JNE1.100 

121 


J2s JB • 2 

ZONEt.lOl 

S2) 


JJ X JB 

EJNEL. 1 02 

i2« 

651 

continue 

Z3NEL.10J 

J24 


JML(Jl.L) * r<lN0(J2L(J2.L).J2UJJ.L)f J2U(J2.L5],J2t(JJfL5n 

23NEL.194 

iSl 


IE (J1.nE.JH) CO 10 66 

2J4EU.135 

ISJ 


Jl* J8N 

ZONEL. 1 06 

554 


J2 » JtJM • 1 

Z3NE4.107 

i5S 


JJ s Jbn ♦ 1 

ZJNEL.I 09 

$5b 


CO TO 651 

ZJN£l. 1 J9 

Jib 

6b 

continue 

ZJNEL.IIO 

Ji6 


iF(LCT,NE.2) CO TO 68 

ZJNEt. 1 1 1 

JbO 


L - LX| 

ZJnEL.112 

J6l 


L2 * L«2 

ZJNEL.ll J 

}6i 


PO 67 J » JNIN.JS 

ZJNEu.lU 

J74 

67 

JNUU.L) * Jf-L(J.L2) 

ZJnEl.I 15 

J76 


CO TO 66 

ZJNEL.llb 

$7T 

6» 

L » L«2 

ZOmEu. 1 1 7 

400 


L2 * t«J 

ZJ > ;l. 1 19 

4 02 

64 

continue 

ZjNit.lH 

402 


DO 74 J » JHIN.J4 

ZJnEl. 1 2u 

4 0 4 


IE ( J.t ). JHN) Cii TO 70 

Z JnEu. 1 21 

4^0 


JP X J ♦ 1 

ZJNi., 1 22 

421 


JTH? I «IN0( JMl( J.L) , JHU( J.L2) , J**L(JP»L) . JNL(JP|L2) ) 

ZJNEU.12J 

4 JO 


IE ( JlMP.fcO. J2LCJ.L) ) CJ TJ 70 

Z3NEU.124 

4 Ji 


JIE (J,L) * 1 

Z3NEL. 1 25 

4 J6 


J2L(J.L) * JTkp 

Z JNE^, 1 26 

4J7 

70 

continue 

ZJNEt.127 

4 42 


IF(LCI.E0.2) CO TO 15 

ZJ<£l.129 

4 44 


nETUnn 

ZJnEl.129 

444 


END 

ZJNEL. 1 JO 
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